Droplet size distribution in homogeneous isotropic turbulence 
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We study the physics of droplet breakup in a statistically stationary homogeneous and isotropic 
turbulent flow by means of high resolution numerical investigations based on the multicomponent 
lattice Boltzmann method. We verified the validity of the criterion proposed by Hinze (1955) for 
droplet breakup and we measured the full probability distribution function (pdf) of droplets radii 
at different Reynolds numbers and for different volume fraction. By means of a Lagrangian tracking 
we could follow individual droplets along their trajectories, define a local Weber number based on 
the velocity gradients and study its cross-correlation with droplet deformation. 



I. INTRODUCTION 



Droplets in turbulent flows occur in variety of industrial processes such as sprays, colloid mills, and mixers. 
When the droplet diameter are smaller than the dissipative scales of turbulence, one can approximate them 
as point particles whose dynamics is governed by the Maxcy-Rilcy-Gatignol equations 1 ^. Several numerical, 
analytical, and experimental studies have studied this regime and found interesting phenomena such as 
clustering of particles at small and large scales '' 4 . When the diameter of the droplet is larger than the 
dissipative scales of turbulence, the point particle approximation is no longer valid. Inertial scale sized 
droplets deform, break, and coalesce under the action of turbulence. The degree of deformation is governed 
by the ratio of the surface tension forces and the intensity of turbulence or the Weber number, We: 

W.m^msM, (1) 
a 

where p^ m ' is the density of the carrier medium fluid, Sujj is the average velocity difference across the 
droplet, the angular bracket indicate spatial averaging, D the droplet diameter, and a the surface tension. 
If inertia forces overwhelm the surface tension We > 1, droplet breaks. Using Kolmogorov theory for velocity 
differences ((Sup) 2 ) ~ e 2 / 3 !? 2 / 3 where e is the energy dissipation rate, Hinze showed, in his 1955 seminal 
work', that the maximum droplet diameter that does not undergo breakup is given by criterion: 

\ " 3/5 
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where the coefficient C = 0.725 is obtained by fitting with experimental data. The above argument does not 
take into account the coagulation of droplets which constitute an important mechanism in dense suspensions. 
Even in dilute limit because of droplet breakup and collision events one expects the droplet dispersion pdf to 
have a finite width peaked around D max . Breaking rate is strongly correlated with the underlying turbulent 
stress coarse grained on a scale of the size of droplet. The latter quantity in fully developed turbulence is 
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strongly non-Gaussian, with maximal deviations from Gussianity observed in the viscous-inertial intermedi- 
ate range. Understanding the stationary droplet distribution due to break-up and coagulation at changing 
Reynolds numbers and droplet volume fraction remains a key unsolved problem. 

Experimental studies by Pacek et al. > focused on turbulent dispersions where the carrier medium and 
the droplet phase have the same density and viscosity and observed droplets size distributions consistent 
with log-normal. Andersson et al. studied both bubble and droplet breakup experimentally and found 
that single bubbles primarily undergoes binary breakups while droplets typically present multiple breakups. 
Risso and Fabrc studied the oscillations and breakup of bubbles in turbulent flow 1 ' while Ravelet et al. w 
focused on breakups of a bubbles in turbulence rising due to gravity. The experimental studies of Eastwood 
et al. focused on the breakup of bubbles in turbulent jets and reported deviations from what predicted by 
the Hinze criterion'. 

To date few numerical studies exist, this is probably due to the additional difficulties implied with the 
need for interface tracking under highly turbulent conditions. Further, to reach good statistical convergence, 
very long and computationally expensive stationary runs are needed. Numerical works by Qian et al. 12 
studied first few breakup events in a high-density contrast bubble breakup and were able to reproduce the 
experimental results of Fabre et al. . Derkscn et al. studied how turbulence history can effect coarsening 
by switching off turbulence in a droplet dispersion. 

The aim of the present paper is to understand the interplay between turbulent fluctuations and surface 
tension in a turbulent emulsion. To this end, we always keep densities and viscosities of the droplet and of 
the medium identical (p^ / p^ — 1 and ^( d )/jy( m ) = 1). We focus on cases with different droplet volume 
fractions 4> = Vd/V m , where Vd and V m design volumes of the droplet or of the carrier phase respectively. 
Furthermore, to avoid any effects of boundaries, we limit to homogeneous and isotropic turbulence. Our 
main results are: (a) we show that a stationary droplet emulsion can be maintained for arbitrary long times, 
with a chosen turbulence intensity and volume fraction, 0; (b) the average droplet diameter is in agreement 
with what predicted by the Hinze criterion at least for small values of 0; (c) droplet pdf show peaks close to 
the average droplet radius and broadens at increasing 0; (d) by means of Lagrangian tracking of individual 
droplets trajectory we show that breakup events are typically associated with peaks in the local energy 
dissipation rate in the droplet neighborhood. 



II. NUMERICAL METHOD 



We model the droplets and the carrier fluid phase by means of a multicomponent Lattice-Boltzmann 
(LBM) algorithm with, by now standard, non-ideal interactions as introduced by Shan-Chen 14-1 ''. This is a 
well established numerical method thus we provide here only its key details. The stirring mechanism, which 
is needed in order to keep the system in a stationary turbulent state, is applied via a large-scale forcing. 

The lattice Boltzmann equations for the Shan-Chen multicomponent D3Q19 model are: 



ft\x + ei ,t + 1) = ft\x) - -[ft\x,t) - ft\p,u)] 

Try 



f\ a \p,u)=p^ Wi 
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Here f"(x, t) are the LBM distribution function at position x and time t for the fluid component a = {d, m} 
(droplet and medium, respectively). The fluid densities and velocities of the individual components are p^ 
and u( a \ Here Wi, are the LBM weights and lattice speeds, i = {0, . . . , 18} 17 ' 18 . The total density of the 
fluid is the sum of the density of the two components p = ^2 a p^ while the total hydrodynamic velocity 
is u = ^2 a «( a '/)( a ' /p. The global effective kinematic viscosity of the fluid is related to the relaxation times 
of its components v = ^2 a c 2 s {t^ — 0.5) 16 , — p^ / p is the concentration, and finally c s = 1/V3 is 
the speed of sound on the lattice for the D3Q19. 

The non-ideal nature of the fluid is introduced by adding an extra force to the LBM equilibrium velocity 
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TABLE I. Runs parameters. For runs N128A - 128D and N512A, p (m) + p (d) = 2.4. For run N512B, p (m) + p (d) = 1.077. 
Here v is the viscosity of the multicomponent fluid, G is the interaction strength and Rx is the Taylor's based 
Reynolds number J1 . The relaxation time of the two fluids r = 0.515 was kept fixed. 
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The non- ideal interaction, as proposed by Shan and Chen, is 15 : 

4c = -Gp {a) (x) P {a ' ) {^ + e l )w l e l (4) 

where {a, a'} = {d, m} and the coupling parameter, G, determines the strength of the microscopic interac- 
tion and effectively sets the value of the surface tension and the diffusivity . Under appropriate conditions 
this force allows the formation of interface between the different fluid components. 

As already anticipated in this study we limit ourself to the case = t^ 71 ' = r which implies v^ d > — 
= v, i.e. identical viscosities for the two phases. The total fluid density + p' m ) can be considered as 
constant over the entire domain, except for the small variations at the droplet interface. This approach to 
model turbulent emulsion has the advantage of being simple and computationally efficient. It was however 
noticed that for extremely long simulation times, as those needed in order to collect firm statistical informa- 
tion, one can observe important diffusion of one fluid component into the other. This physical phenomenon 
can be artificially large in the LBM simulations due to the relatively big interface width between the two 
fluid components. We have however shown in a recent publication-'" how this effect can be controlled by 
means of an effective and computationally efficient manner. By means of this improvement, the presented 
algorithm allows us to investigate turbulent emulsions for arbitrarily long simulations times. 

In order to stir turbulence we apply forcing at each position and at each time step modulated by means of 
a sum of sine waves with small wavenumbers. To ensure an homogeneous and isotropic stirring, the phases 
of the sine waves are evolved in time by means of a stochastic process. The forcing is divcrgcncclcss and its 
expression for the generic i-th component is: 

where i, j = {1, 2, 3}, F^ Q ' is the external force at x and time t. A controls the forcing amplitude, kj are the 

wave- vector components and the sum is limited to k 2 = k\ + k\ + fc| < 2. The phases are evolved in 
time according to independent Ornstein-Uhlenbeck processes with the same relaxation times T = u rms /N 
(N is the linear domain size and u rms was taken equal to 0.1 to almost match with the typical values for 
the large scale velocity). 
III. THE STATIONARY STATE 



In presence of turbulence droplets continuously undergo breakup and collision processes, whose rate 
depends on both volume fraction and Reynolds number. The plot in Fig. 1 shows a typical time evolution of 
the droplet number for volume fractions (f) equal to 0.5%, 5% and 10% at Re\ = 15 (runs N128B — D, Table I). 
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5 t/x eddy 15 20 25 

FIG. 1. Number of droplets vs. time for different values of the volume fraction, <f>. An increase in the volume 
fraction leads to an increase in the total droplets number and more coagulation processes occur. The time is made 
dimensionless by the large scale eddy turnover time, r e dd y = N/u rma . Top panel shows snapshots of the droplet 
dispersion at the intermediate t/r e ddy = 0.33 and steady state t/r e ddy = 12.5, 25 configurations. 



During the early simulation stages, when turbulence is still developing, the droplet deforms and then starts 
to breaks up into smaller droplets. At later times, the droplet count steadily oscillates around a mean value 
because of a continuos competition between breakups and coalescence events. The stationary probability 
distribution function of droplet sizes, with its corresponding cumulative distribution, is shown in Fig. 2 for 
different volume fractions <j> = 0.07%, <fi = 0.5% and <f> = 5% with Re\ = 15 and a = 1.6 • 10~ 3 (runs 
N128A- C, Table I). 

In the case 4> = 0.07% only very few breakups are observed and the droplet PDF closely resemble a delta 
function. We observe a mild increase in the average droplet radius at increasing the volume fraction of the 
droplet phase; furthermore, because of enhanced coalescence, we observe that the maximum droplet radius 
increases with increasing the volume fraction (see Fig. 2 (bottom)). In the same figure we also superpose the 
empirical pdf with a log-normal distribution. 
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As one can see, the two fits are in good agreement with the data, within statistical errors, at large and 
intermediate radiuses. At very small droplet radius one can observe deviations between measured data and 
the fitting function. On this point we are not able to draw conclusive statement and we tend to think that 
deviations may be unphysical and due to the limitations of the diffused interface models to capture droplets 
whose size is comparable to the interface width. 
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FIG. 2. (top-panel) Log-log plot of the probability distribution function (pdf) of droplet radius for different values 
of the volume fraction. Pdf are scaled to have maximum value unity, (bottom-panel) Cumulative probability 
distribution function for different volume fractions. In panels, dashed lines represent log-normal fits to distributions 
with = 0.5% and <j> = 5%. 
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FIG. 3. Probability distribution function of the droplet radius from two simulations with almost the same volume 
fraction but with different Re\. The inset shows the pdf normalized by the average droplet radius. 



To understand the effect of the Reynolds number on the stationary droplet size distribution, in Fig. 3 
we compare the droplet radius PDFs with comparable volume fractions but for different Reynolds number 
(j) = 0.5%, Re x = 15 and <j) = 0.3%, Re x = 30 (corresponding to runs N128B and N512A, Table I). Our result 
support the validity of the Hinzc criterion which predicts that, for the same material parameters, the average 
droplet diameter should only depend inversely on the energy dissipation rate [Eq. (2)]. By rescaling the 
droplet radius with their average values, the two PDF collapse (inset of Fig. 3). Thus, we can exclude strong 
dependency of the pdf shape, at fixed volume fraction and material properties, on the Reynolds number; at 
least this is the case for the two Reynolds number that we investigated. 



IV. HINZE CRITERION 



The Hinzc criterion, as in Eq. (2), provides the estimate for the maximum droplet diameter, D max , that 
should not undergo breakup at a fixed turbulence intensity. In presence of turbulence however breakup and 
coalescence always occurs and the correct indicator for the droplet size is the average droplet radius. We 
replace therefore D max by the average droplet diameter (D), where with the angular bracket we indicate 
averaging over a steady distribution of droplet diameters. From Fig. 4 one can see that our data support 
the validity of Hinzc criterion in the case of small volume fractions, while we one can clearly detect a non 
trivial dependency on the volume fraction for larger values of <f> (keeping the turbulent intensity and surface 
tension constant). 

As the droplet is transported by the turbulent flow, it visits regions of where velocity gradients are much 
larger or smaller than the average value. Thus, a better characterization of the droplet breakup required the 
tracking of individual droplets along their trajectory along with monitoring their deformation and the local 
flow conditions in the immediate neighboroud, i.e. the local Weber number, We s . Here wc define the local 
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FIG. 4. Droplet diameter versus the energy dissipation rate to test the Hinze criterion for droplet breakup '. The 
different symbols correspond to runs that are distinguished by their volume fraction <f> and surface tension a (see 
Table I). The solid black line is the prediction of Hinze [Eq. (2)]. Note how at larger volume fractions deviations 
from the Hinze criterion are observed, these are probably due to enhancement of coagulation events. The error-bars 
indicate the uncertainty in the estimation of the droplet interface position. 



Weber number as the difference between the maximum velocity difference between a point on the droplet 
surface and the velocity of the droplet center of mass We s = \u maXl s — Ucm\- The droplet deformation is 
instead characterized by a dimensionless parameter S* = (S/Sq — 1) where S is the surface area of the droplet 
and So = (4 7 r) 1 / 3 (3y d ) 2 / 3 is the surface area of the equivalent spherical droplet (i.e. a spherical droplet with 
the same volume of the deformed one). With this definition S* = corresponds to an (undeformed, i.e. 
spherical) droplet, whereas S* > indicate a stretched droplet. The plot in Fig. 5 shows the time evolution 
of the surface deformation, S*, of the local Weber number We s , and times at which breakup event occur, 
during a part of the droplet trajectory. Data in Fig. 5 refers to the run with the smallest volume fraction, 
<f> = 0.07, and with Re\ = 15. For this run, due to its low droplet count, it was technically easier to track 
breakup events (run N128A, Table I). As one can see, the droplet deformation, S*, correlates strongly with 
the local Weber number, We s . Furthermore, just before the droplet breaks (stars in the figure) one can 
identify a sharp increase in both surface deformation S* and We s . We can further quantify the correlation 
by computing the joint pdf of S* and We s over the full time-series for run N128A (see Fig. 6). From this 
joint PDF we find that S* ~ Q.QQQWe. In Rcf. ~ it was reported a similar correlation but for the rather 
different case of a vapor-liquid dispersion with a density contrast between the two phases. 
V. CONCLUSIONS 



We have studied the droplet dispersion in a turbulent emulsion from very small to larger volume fractions. 
We focused on the role of turbulence fluctuations on the droplet deformation and breakup and to this end we 
limited to density and viscosity matched emulsions. In this way we could better highlight the new physics 
associated with the interplay between surface tension and turbulence fluctuations. Our numerical study, 
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FIG. 5. Correlation between the Weber number and the non-dimensional droplet extension S* = (S/So — 1). The 
instantaneous Weber number is computed by using the maximum value of the energy dissipation at the droplet 
surface We B . The plot represent the evolution of a typical droplet trajectory from the run with smallest volume 
fraction <j> = 0.07% and Re\ = 15 (N128A). Stars indicate the breakup events that we observed during this particular 



which is based on the lattice Boltzmann method, shows that by using a proper large-scale isotropic stirring 
mechanism 20 , it is possible to attain a stationary steady state which allows to study the physics of droplets 
distribution and evolution in turbulence. The pdf of droplets radii follow a log-normal distribution and the 
Hinze criterion is well satisfied at small volume fraction, while we observe a departure from its prediction 
at higher volume fractions. We performed a Lagrangian tracking of individual droplets and we showed that 
the local Weber number is a suitable prognostic indicator for droplet breakup and correlates strongly with 
droplet deformation. 
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